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I. INTRODUCTION 



New sources of short radiationpulses, in particu- 
lar high-harmonic generation 0, Q and free-electron 
lasers Q, promise to enable a new generation of 
pump/probe experiments on molecules in which the cen- 
tral frequencies of both pulses are in the ultraviolet or 
X-ray. The short time scales that have recently become 
practical for these measurements extend to subfemtosec- 
ond pulses with delays between them on the order of 
femtoseconds or tens of femtoseconds. This generation 
of experiments can also involve probe pulses which ion- 
ize or dissociate the target molecule and thus add the 
dimension of the time-resolved measurement of the elec- 
tron and molecular fragment energy distributions. To 
accurately interpret and describe the results of these ex- 
periments, ab initio methods must be developed to treat 
highly electronically excited and strongly nonadiabatic 
molecular dynamics. Coincidence experiments for in- 
stance, demand the capability to describe multiple ion- 
ization and dissociation, and nonlinear effects [H, Q may 
entail the excitation of multiple electrons. In general, 
most of these phenomena currently remain beyond the 
reach of accurate ab initio theoretical descriptions. 

However, significant headway has been made. Ap- 
proaches employing classical trajectories on coupled 
Born-Oppenheimer surfaces @, [1| are well suited for dy- 
namics on low lying excited bound electronic states, and 
have been applied to molecules as large as DNA bases @ . 
With an approximate treatment of the coupling to the 
ionization continuum, such trajectory methods may de- 
scribe time-resolved photoelectron signals in pump-probe 
experiments 0, Q ; such a treatment also permitted the 
study of the quantum nuclear dynamics of Auger de- 
cay M- 

While these methods have already shown great 
progress in describing a range of non-Born-Oppcnhcimer 
and excited-state effects, their utility is greatest for situ- 



ations in which ionization may be treated approximately 
and in which the interacting electronic states may be ex- 
plicitly identified. Several other approaches avoid the use 
of Born-Oppenheimer states altogether and show promise 
for treating highly electronically excited or nonadiabatic 
electronic and nuclear dynamics [ill . [T3 |. In a similar 
context, ionization has been included in a variational 
treatment that explicitly includes electron-nuclear cor- 
relation |l3| and has also been treated with coupled elec- 
tronic and semiclassical nuclear wave packets [14j . 

The time-dependent multi-configuration Hartree Fock 
(MCTDHF) approach would seem to be a natural and 
widely applicable starting point to the electronic part of 
this problem, because it is capable in principle of ex- 
actly describing the dynamics of many-electron motion. 
There is a considerable literature on this subject already 
including analysis of the formulation of MCTDHF and 
its application to small or model systems (l5rl23| . and 
also attacks on the combination of electronic and nu- 
clear motion j24|. On the basis of this literature, the 
fundamental idea can be said to be well established that 
a time-dependent linear combination of determinants of 
time dependent orbitals should be flexible enough to de- 
scribe the electronic response a molecule to short intense 
pulses in any part of the spectrum. 

Similar ideas in a different context have under- 
pinned the development of the multiconfi gur ation time- 
dependent Hartree (MCTDH) method [25M30I ] . which has 
had considerable success in treating problems of nuclear 
dynamics, including vibronic coupling and reactive scat- 
tering. However, by comparison, the MCTDHF method 
for electrons has still not delivered its full potential, par- 
ticularly in the presence of ionizing fields or including 
nuclear motion. The reason appears to lie in several seri- 
ous technical barriers to its implementation and general 
application: 

• Electronic and nuclear motion are strongly correlated; 
the cusps in the electronic wave function at the posi- 
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tions of the nuclei must be accurately represented for 
all geometries, and the basis set error in the electronic 
part of the calculation must not depend strongly on 
nuclear coordinates. 

• The evaluation of the two-electron integrals over the 
time-dependent orbitals must be numerically efficient, 
otherwise it will dominate the computational time re- 
quired. 

• The ionization continuum must be properly treated 
within the MCTDHF description, if it is relevant to 
the problem at hand. 

• The nonlinear, unitary, stiff differential equation in- 
volving orbitals and configuration coefficients must ef- 
ficiently numerically integrated. 

Here we address all these difficulties for the case of di- 
atomic molecules. The key to overcoming the first three is 
the use of the Finite Element Method and Discrete Vari- 
able Representation (FEM-DVR) in prolate spheroidal 
coordinates. The electronic basis is a set of piecewise 
interpolating polynomials with cusps at the nuclei, para- 
metrically dependent upon the bond distance. Such an 
atom-centered, parametrized basis has been used be- 
fore [13j], but the choice of prolate spheroidal DVR has 
several important advantages. 

• Orbitals defined in these coordinates are orthogonal 
for all values of the internuclear distance R, and there- 
fore a single set of orbitals may be used for all nuclear 
geometries, radically reducing numerical effort. The 
prolate spheroidal DVR basis allows for a sparse rep- 
resentation of the primitive two-electron integrals and 
their rapid contraction into orbital matrix elements, 
and it leads to rapid convergence of both bound and 
continuum wave functions, as we have found in previ- 
ous fixed-nucei calculations on single and double pho- 
toionization of H2 (3ll - [33j . 

• The basis enables rigorous inclusion of the ionization 
continuum via Exterior Complex Scaling (ECS) [3^ |. 
The implementation of the ECS formalism with the 
FEM-DVR approach is well established as a formally 
sound and computationally efficient treatment of both 
photoionization and electron- impact ionization [35j . 
because it imposes correct outgoing scattering bound- 
ary conditions for both single and double ionization 

• Finally, inclusion of nuclear motion in the case 
of diatomics, without the necessity of the Born- 
Oppenheimer approximation, is also greatly simplified 
by the prolate spheroidal FEM-DVR. We employ a 
good approximation to the exact nonrelativistic molec- 
ular Hamiltonian including the interaction with the ra- 
diation field, omitting only the Coriolis coupling and 
mass polarization terms. 



Thus, the prolate spheroidal DVR basis is ideally 
suited for the study time dependent excited state dynam- 
ics of diatomic molecules. Its matrix elements appear in 
a nonlinear differential equation of potentially large di- 
mension, however, and the integration of this equation 
is a formidable challenge in its own right. Several differ- 
ent methods of integrating the MCTDHF equations have 
been described, and we have found an efficient and stable 
generalization of the method of Ref. [27j . 

The choice of a single set of electronic orbitals with 
parametric dependence on R that is used in this ap- 
proach might appear to be, from the traditional perspec- 
tive of the Born Oppenheimer approximation, an unnat- 
ural starting point, and raises questions regarding the 
convergence of the wave function for coupled electronic 
and nuclear motion with respect to the numbers of con- 
figurations included, which we must address. However, 
this choice also offers important advantages, simplifying 
the working equations in a critical way and allowing the 
slow, nuclear degree of freedom to be distributed across 
supercomputer processors. 

The outline of this paper is as follows. In Sec. UH wc 
first review the working formalism for MCTDHF with- 
out nuclear motion. In Sec. Mil wc describe the formula- 
tion of the diatomic problem in prolate spheroidal coor- 
dinates, and the form of the Hamiltonian appropriate to 
those coordinates when the underlying electronic basis 
has a parametric dependence on internuclear distance. 
Sees. IIVI and [V] describe the inclusion of nuclear mo- 
tion using a DVR basis in the nuclear coordinate R in 
combination with the MCTDHF treatment of electronic 
motion. We then discuss the application of the ECS 
transformation to treat ionization within the MCTDHF 
framework in Sec. IVII The remaining sections discuss 
computational details and some preliminary results for 
both bound and continuum electronic motion. We use 
atomic units throughout. 



II. MCTDHF FORMALISM FOR FIXED 
NUCLEI 

The MCTDHF working equations have been formu- 
lated previously fl5l - [l7l . list [2a l23j , and we will give here 
only a brief description of the working equations in or- 
der to establish the starting point for the inclusion of 
nuclear motion using this approach, and to indicate how 
exterior complex scaling of the electronic coordinates is 
implemented in this context to treat ionization. The 
MCTDHF approach begins with an expansion of the elec- 
tronic wave function in antisymmctrized products (deter- 
minants) of time-dependent spin orbitals 



|*(t)>= $> a (i)|n a (*)> 



(1) 



m 
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which each antisymmetrized product of N spin orbitals 
specified by the vector ft a and is defined by 

\n a (t)) = (\<f> nal (t)) x ...\c/> naN (t))) . (2) 
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We use spin restricted orbitals which are product func- 
tions of the space and spin coordinates % = {r i; S,} of a 
single electron, 



(q\(f>n(t)) = <j> an (r) TOS „ (5) 



(3) 



where Q is a spinor, Q = a or Q 
second quantization we can write 



Alternatively, in 
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in terms the creation operator, a', for spin orbitals. The 
space part of each spin orbital is expanded in a set of 
time- independent DVR basis functions, fj, which we will 
specify in Sec. IIIII below. 



(5) 



where the projector P is the matrix representation of 
the projection operator, P = ^ Q \(f> a (t))(4> a (t)\, onto the 
space spanned by the orbitals at time i, 



(9) 



so that 1 — P projects this equation on to the space or- 
thogonal to that spanned by the orbitals. Our conven- 
tion is that boldface symbols are matrices in either the 
orbital (c) or configuration (A) basis, and Ov stands for 
matrix multiplication. In Eq. ([8]), p Q7 is the reduced 
one-electron density matrix for the wave function in Eq. 

©■ 
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(10) 



associated with a timc- 
c a (t); the vector of all c a 



so each spatial orbital is 
dependent coefficient vector 
is denoted c. 

The MCTDHF equations are based on the application 
of the Dirac-Frenkel variational principle for the time- 
dependent Schrodingcr equation to the trial function in 
Eq. (J), 
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where we include the constraint that the orbitals remain 
orthonormal, (4> a \4>p) — 5 a p = 0, along with the corre- 
sponding Lagrange multipliers \ a p. The variations in 
Eq. ([6]) are variations in the coefficients A^ and c a . 

In this work we employ a full configuration interac- 
tion (CI) representation of the wave function in Eq. (fT]). 
although further application of the method could entail 
restricted CI wave functions such as the "complete active 
space" approach that are standard in modern quantum 
chemistry. An important point is that in the full CI case 
addressed here the solution of Eq.© yields A Q( g = 0. 
The wave function is invariant with respect to rotations 
among the orbitals, which may be compensated for by ro- 
tations among the A-coefficients. The solution of Eq.© 
is therefor not uniquely defined and unique time propa- 
gation requires an additional constraint besides orthog- 
onality of the orbitals. The simplest constraint is to set 
the time-derivative matrix g in the orbital basis to zero, 
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For the orbitals one obtains the equation of motion 



^c a = ( i-p)E 



L 



C0, (8) 



and is the matrix of one-body operators in 

the Hamiltonian with respect to the underlying time- 
independent basis in Eq. ([S]). All quantities in Eq. 
© are time-dependent except for the identity and h^ 1 ) 
matrices (in the absence of an external time-dependent 
field). The reduced two-electron operator W is de- 
fined [TU, [lH in terms of the reduced two-particle density 
matrix, r 7sa ; and the electron repulsion, W s ; , expressed 
as its matrix representation in ri in the underlying time- 
independent basis, 



W7/3 = ^E r T-«' w «»(*) 



w Bl (t) = / rM,t) 



(ii) 



R 



In - r 2 \ 



-Mr 2 ,t)dr 2 . (12) 



In Eq. p^]) the electron repulsion operator appears as 
R/\fi — r 2 \ because its matrix elements then have no R 
dependence in the underlying prolate spheroidal DVR, 
and therefore also in the orbital basis. That fact signif- 
icantly simplifies the implementation of the MCTDHF 
equations in these coordinates. 

The time derivative of the orbital coefficients being 
specified, one then obtains the equations of motion for 
the A-coefficients 
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(ri a \H\ri a > 



(13) 



where H is the matrix of the Hamiltonian in the config- 
uration basis, consistent with Eq.©. 

For Hermitian Hamiltonians, the MCTDHF equations 
conserve the norm and the expectation value of the en- 
ergy [26( . We next turn to the specification of the under- 



lying basis in Eq. © in prolate spheroidal coordinates and 
the implementation of ECS in that basis for the treat- 
ment of ionization during the propagation of Eqs. (|13[) 
and ©. 
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III. FIXED-NUCLEI HAMILTONIAN AND 
WAVE FUNCTION 

A. Hamiltonian 

Prolate spheroidal coordinates were used in early quan- 
tum chemical calculations on diatomic molecules, be- 
cause analytic basis sets, such as Slater-type orbitals, in 
those coordinates could exactly satisfy the cusp condi- 
tions on the two nuclei, and because their scaling proper- 
ties with internuclear distance offered additional compu- 
tational advantages [U |42| . We use them here for some 
of the same reasons, but also because the implementation 
of the FEM-DVR approach in these coordinates dramat- 
ically simplifies the calculation of the two-electron inte- 
grals. 

If the nuclei of our diatomic molecule are at Ra and 
Rb, we can define prolate spheroidal coordinates (£, ij, ip) 
for each electron in the usual way by rotating a two- 
dimensional elliptical coordinate system (£, rj) about the 
focal axis of the ellipse, 



'/ 



\r-R A \ 



Rb\ 



R 

R A \-\r-R B \ 
R 



(1 < £ < oo) 
(-1 < i? < 1) 



(14) 



and the remaining coordinate, ip, (0 < ip < 2ir) is the 
azimuthal angle. The one-body operators in our Hamil- 
tonian in these coordinates are specified by the Laplacian, 



R*{e-if) 
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and the electron-nuclear attraction 
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(15) 



(16) 



(17) 



\r-R A \ \r-R B \ - R(e 

while the one-electron volume element is 

dV = {R/2f{i 2 -f] 2 )d£, dr] dip . 

For an N-electron problem then, a factor of i? 3 appears 
in the volume element in Eq. (|17[) for integration over the 
coordinates of each electron. We eliminate that factor 
in a fixed-nuclei calculation by solving for R 3N / 2 times 
the electronic wave function. In Sec|V]when we consider 
nuclear motion we will solve for i^ 3JV / 2 + 1 times the total 
wave function to further simplify the form of the nuclear 
kinetic energy operator. 

B. Wavefunction 

We make use of an FEM-DVR in both £ and rj for 
each electron, and our primitive basis functions are prod- 
ucts of those DVR functions with a factor describing the 



4> motion with a particular angular momentum projec- 
tion, M, along the molecular axis. As we have noted in 
previous studies [H, on H 2 , the specification of the 
FEM-DVR depends on whether M is even or odd, in or- 
der to properly represent the analytic dependence of the 
wave function as £ or 77 approach the singularity at ±1. 
Since we are constructing a DVR in each finite element, 
we use Gauss-Radau quadrature on in the element in £ 
beginning at £ = 1, and Gauss-Lobatto quadrature for 
the others. We use Gauss-Legendre quadrature to define 
the DVR in rj. So defining the basic DVR interpolating 
functions in terms of the quadrature points and weights 
by (with x — r\ or £) 



1 N 



■ , . %i Xj 



lx 2 - 1 

xT^T 



even M. 



odd M 



(18) 

We define the one-electron primitive FEM-DVR func- 
tions according to 



f M = 

J ia 



a iM<p 



(19) 



Note that if we use the underlying quadrature to cal- 
culate the overlaps, because for fixed nuclei we solve for 
R? N I 2 times the wave function, these functions are or- 
thonormal with respect to the volume element |(£ 2 — 
rj 2 ) d£ dr] dip. 

We have discussed the simplifications of the integrals 
of l/r 12 in a similar basis before (33[, in which we used 
spherical harmonics in the variables rj and if instead of 
the DVR in 77 we are using here. In this basis the sim- 
plifications are even more powerful. By making use of 
the Neumann expansion of I/V12 in prolate spheroidal 
coordinates, and solving Poisson's equation in £ in the 
FEM-DVR basis followed by using the Gauss quadrature 
in rj we arrive ultimately at expressions for the two elec- 
tron integrals of the general form 



1 



n.2 



(20) 



Sis'Sjj'S aM 'Sb : b'Fij !a ,b 



Where in addition we have the requirement that 
m = M 3 - Mi = M 2 — Mi. We give the explicit for- 
mula in Appendix [EJ We can see therefore that the two- 
electron integrals in the FEM-DVR basis are diagonal 
in the indices corresponding to both £ and ry. This is 
an immense simplification over the use of analytic basis 
functions, like Gaussians, and dramatically reduces the 
effort in transforming the two electron integrals from the 
FEM-DVR basis to the time-dependent orbital basis, and 
thereby simplifies and speeds up the both the construc- 
tion of the two-electron portion of reduced Hamiltonian 
in Eq. (fT2"l) and its operation in Eq.©. 
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IV. HAMILTONIAN AND WAVE FUNCTION 
FOR NUCLEAR MOTION 

A. Hamiltonian 

When we include nuclear motion, we must account 
for the fact that this FEM-DVR basis is i?-dependent 
through the dependence of the electronic coordinates on 
internuclear distance. The DVR quadrature points that 
define the basis in Eq. (fT9]) move with R. As has been dis- 
cussed before - for example by Esry and Sadgepour [43| - 
in a basis of functions of prolate spheroidal coordinates, 
the derivatives with respect to R in the Hamiltonian must 
be calculated holding the cartesian coordinates, x,y and 
z of the electrons fixed instead of holding £, 77 and <f> fixed. 
The relation between those derivatcs is 



d_ 

dR 



N 1 

Y- 



R 



Y 



(21) 



Y 



where the sum is over electrons, and the operator Y is 
1 ( d 

(22) 

'drjj 



with a = (Ma — Mb)/ {Ma + Mb) being the mass asym- 
metry parameter. We must use these relations when con- 
structing the second derivative term, (d 2 /dR 2 ) xyz , in 
the nuclear kinetic energy. When including nuclear mo- 
tion, we calculate R 3N / 2+1 times the wave function, and 
so the radial part of the nuclear kinetic energy operator 
then becomes 



K R 




£v<t> 



1 3N\ 



2 J \R \dRj 



iv<t> 



(23) 



For more than one electron, we do not employ an exact 
treatment, which might be accomplished in polysphcri- 
cal coordinates 0, EH for example, but instead employ 
a straightforward adaptation of the one-electron Hamil- 
tonian [43l ] that omits several minor terms unimportant 
for a host of nonadiabatic dynamics relevant to attosec- 
ond physics. In particular, our electron position vectors 
are represented in prolate spheroidal coordinates relative 
to the center of mass of the nuclei. We therefore omit 
terms relating of mass polarization of the electrons, i.e., 
for N electrons, the deviation of the center of mass of 
the nuclei from the center of mass of the nuclei and any 
N-l subset of the N electrons. In the chosen coordinate 
system, such terms would be represented as two-electron 
derivative operators. In contrast, an exact polyspheri- 
cal treatment such as that using heliocentric Radau co- 
ordinates (46j . in prolate spheroidal form, would entail 



nonseparablc corrections to the electron-nuclear poten- 
tial which, without further approximation, would be in- 
tractable in the present framework. 

We additionally omit Coriolis coupling and write the 
Hamiltonian for rotational quantum number J as 



H = Kr 



1 



2 m R2 



J(J + 1) - 2 J 2 + V 



1 

2^ 



(24) 

The reduced masses are defined as /ir = mAVtiB / + 
ttib) where A and B are the two nuclei, and (in atomic 
units) fi e for an N-electron system is fi e = (niA + ms + 
N — l)/(rriA + tub + N). J z is the projection of the 
electronic angular momentum on the internuclear axis 
and equals the sum over the l z eigenvalues of the in- 
dividual electrons, X)i!=i ^i- The operator P is the 
square of the electronic angular momentum operator. 
The potential, V , includes the electron-electron repul- 
sion, electron-nuclear attraction, internuclear repulsion, 
and time-dependent dipole interaction term, if an ex- 
ternal field is being applied. At present, we also omit 
the two-electron terms in I 2 (proportional to Ulj and in 
Yi) 2 (proportional to YiYj). These terms are in gen- 
eral similar in magnitude to the mass polarization terms 
that we have already omitted, and also not relevant to 
the processes of immediate interest. Usual nonadiabatic 
effects such as curve crossing transitions are driven by 
the cross term in Eq. (|23p involving products of electronic 
and vibrational momenta, as opposed to the terms that 
we have omitted containing terms quadratic in the elec- 
tronic momenta. They and the Coriolis coupling may be 
included in future applications, and their omission here 
accelerates the numerical implementation. 



B. Wave function 

To include nuclear motion and electronic motion si- 
multaneously, and also avoid the Born-Oppcnheimer Ap- 
proximation, we begin with a trial function in which we 
use the same configurations specified in Eq.([T|), express- 
ing the time-dependent orbitals in the prolate spheroidal 
FEM-DVR, and taking advantage of the implicit depen- 
dence on the internuclear distance R of those orbitals, 



<fl|*(t)> = 5^,«(*)|tf«(*;i2)>x«(fl). 



(25) 



In Eq.fl2S} the function Xk(R) is an ordinary FEM-DVR 
basis function based on Gauss-Lobatto quadrature within 
finite elements in R, and is labeled by the grid point 
R K at which it is nonzero. The coefficients, A a ^ K {t), of 
the configurations depend explicitly on the index k as 
well as the configuration index a. So this representation 
of the wave function uses MCTDHF for the electrons 
while using a full primitive basis DVR representation of 
nuclear motion. The configurations have a parametric 
dependence on R, which we emphasize with the notation 
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\n a (t;R)), and that dependence comes entirely through 
the R dependence of the prolate coordinates, and thus of 
the FEM-DVR grid for the electrons, 

(q\<l> n {t;R)) = (t> an (r(R),t)n mSn (E) (26) 

where r(R) = (£ i (R),ri(R), ip) are the prolate spheroidal 
coordinates. 

Using this trial function in the variational principle 
in Eq.© means that a single set of orbitals is used to 
describe the electronic motion for all R, and thus that 
the coefficients cf (t) 



(27) 



do not depend on R. This fact simplifies the resulting 
MCTDHF equations in a fundamental way, and is key to 
the practically of this approach. The ansatz in Eq. ([25|) 
is largely motivated by this fact, but it is clear nonethe- 
less that it is capable in principle of representing the ex- 
act wave function if a sufficient number of orbitals is in- 
cluded. Our numerical tests below will test the efficiency 
of this expansion. The coefficients of the configurations 
depend explicitly on the R index, and thus are capable of 
weighting the configurations constructed from those or- 
bitals differently at different internuclear distances, and 
thus different orbitals can contribute differently at differ- 
ent values of R. The nuclear cusps of the wave function 
are represented accurately at all values of R in this ap- 
proach, and the orbitals are automatically orthogonal for 
all internuclear distances. 



V. MCTDHF WORKING EQUATIONS 
INCLUDING NUCLEAR MOTION 

Using the following notation for combining parts of the 
Hamiltonian in Eq ([24|) : 



T 



E 



1 d 2 



V?- 



D' 



— U 

^PR 



E 



Y 



2dR 2 

(28) 

the Hamiltonian including a radiation field employed here 
for R^ N / 2 + l times the wave function, omitting Coriolis 
and two-electron terms in I 2 and Yi) 2 , may be com- 
pactly expressed as 



R R hr V 



T R + jjRjjel 



H{t)=H +R£{i)fi (length) 
1 A{t) 



R c 



-fi (velocity) 



(29) 



(30) 



where £(t) and A(t) are the electric field and vector po- 
tential, and (jl is a coordinate or derivative operator in 
the electronic prolate spheroidal coordinates. The oper- 
ators D R and D el are antihermitian, first order differen- 
tial operators, and the potential is a separable product of 
1/R times a potential for the nuclear repulsion plus two 
electron repulsion that is a function of only the prolate 
spheroidal coordinates: 

V = (Z X Z 2 +«12(£lj 7?l,<£l,£2, V2,(P2)) (31) 

with V12 = R/r\2- 

The working equations for MCTDHF with nuclear mo- 
tion are similar to the Born-Oppcnhcimcr version. The 
inverse of the reduced single-particle (electronic) reduced 
density matrix still appears, and that matrix is a sum 
over nuclear grid points, 



Pa/3 



^ , PnaK(3 7 



(32) 



where 



PH.aT/3 = ^2 A l,K A b,T{n a \alims a Pms\nb) , (33) 
a b ms 

and the reduced two-electron matrix W 7 ^ is defined sim- 
ilarly. We also have the reduced two-electron density ma- 
trix defined as 

q a' k /3 (3' = 

E K,K A bAna\al m s a l'm s > a 0rnsap> ms >\n b ) . 
a b ms ms' 

(34) 

By defining reduced matrices for 1/R, 1/R 2 , and R 
(which appears the length gauge dipole operator), and 
a reduced derivative operator, 

M< S = E' 9 ««^ i? « Q = 1.-1,-2 , 



^ t r« ol ol ! k, j3 f3' j-y j 

K 

? _ \ ^ n n R 

a/3 — / Pk a r (3*-* kt ; 



(35) 



along with matrix elements of the differential operators 
in Eq. (|2"5|) . we arrive at the MCTDHF equations for the 
orbital coefficients, 



d 



^ = E( 1 - p )^E 



7 



7>{ — 2)rpel 

^7 L 



(36) 



The equation for the coefficients of the configurations still 
has the form ij^A — HA. We provide expressions for 
the various matrix elements of the one-electron operators 
appearing in this section in Appendix [X] 
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VI. EXTERIOR COMPLEX SCALING AND 
THE TREATMENT OF IONZATION 

In the application of the MCTDHF approach to 
molecules subject to short UV pulses, there is in gen- 
eral some component of the wave function that is ion- 
ized. As the wave function is propagated the outgoing 
electron flux will inevitably reach the end of the FEM- 
DVR grid in £ and reflect back from it. In a number 
of other time-dependent applications of grid methods it 
has been shown that the ECS transformation is capable 
of perfectly extinguishing those reflections, both in the 
absence of an external field (4lJ and in the presence of an 
time- varying field [H, [49| . 

In a time-dependent calculation the ECS transforma- 
tion enforces outgoing wave boundary conditions for the 
ionized part of the wave function, even at very long 
times. In the application of the ECS method in pro- 
late spheroidal coordinates the electronic coordinates are 
scaled only beyond a radius £o by a complex phase factor 
according to £ -> £ + (£ - £,o)e ie , where < 6 < ir/2 
is an angle on which the results do not depend formally. 
The value of £o is chosen large enough that the physi- 
cal quantities of interest can be calculated from the wave 
function for all electronic coordinates satisfying £ < £o- 
This is a formally exact procedure, and in a converged 
calculation does not alter the wave function in the inner 
region from its exact value. In fact, we may calculate ac- 
curate bound state energies even with £o = 1, as shown 
in Table HI However, to extract ionization information, 
we choose a larger £o and perform analysis inward of that 
value, on the real £ axis. 

The analytic continuation of the Hamiltonian under 
the ECS transformation leads to a complex symmetric 
matrix representation when the basis functions are real 
at real values of the coordinates, and the DVR implemen- 
tation of ECS, detailed previously [HI, [5(| also leads to 
complex symmetric matrices. So for example, the matrix 
representation of the one-body operators in the Hamilto- 
nian, hW in the present FEM-DVR is complex symmet- 
ric, and the DVR basis functions, as defined in Eq. (fT8|) 
are orthonormal when their overlap is quadratured along 
the ECS contour 

Once the MCTDHF orbitals are expanded in terms of 
the orthonormal FEM-DVR basis in Ed-lfT^)). they are 
represented by the coefficient vectors c a in Eq. (|27|l . We 
may define the inner product of a pair of those orbitals 
in terms of the expansion coefficients in the usual, Her- 
mitian way, 

(0a 1 0b) = c a f • c b (37) 

making use of the orthonormality of the DVR functions 
on the ECS contour. We can then arrange the coefficients 
defining the MCTDHF orbitals as the matrix Ci yCC , and 
use this inner product when transforming the operators 
from the FEM-DVR basis to the orbital basis. So for 
example we would transform the ECS-scalcd one-body 



Hamiltonian according to 

h« = Cth« s C (38) 

where f denotes the Hermitian conjugate of the matrix. 
Because it is unitary (in the limit that there are the 
same number of orbitals as FEM-DVR basis functions) , 
this transformation does not change the spectrum of the 
ECS-scaled matrix representation of h^. Doing every 
transformation to the orbital basis that is involved in 
constructing the matrices in the working equations, Eqs. 
t[13]> and ([8]). in this way preserves the analytic properties 
of the ECS solutions. 

The implementation of ECS in this manner in the 
MCTDHF equations has another, very important advan- 
tage. As the solutions are propagated forward in time, 
the constraint that the orbitals remain orthogonal, and 
the constraint on their variations in time imposed by 
Eq.([7]), are then imposed with the usual sense of the 
Hermitian inner product in Eq. (|37l) . This procedure is 
essential to the numerical robustness of the MCTDHF 
method, because if an inner product without complex 
conjugation were used instead, as it is frequently in the 
complex scaling literature, an orbital could in principle 
have zero overlap with itself. While there may be no for- 
mal reason to choose one implementation over the other, 
there is therefore a compelling numerical reason to chose 
the Hermitian inner product. This choice does result in 
matrices H and W that appear in the working equa- 
tions, Eqs. (|13[) and (JSJ for example, with no symmetry, 
but nonetheless the overall properties of the solutions, in 
particular their outgoing wave character, under ECS are 
correctly reproduced. 

VII. NUMERICAL INTEGRATION 

The integration of the coupled nonlinear differential 
equations of Eqs. tfTB"]) and © has been the source of sev- 
eral theoretical and numerical studies over the past years. 
Splitting of the orbital equation by separating the one- 
body, stiff kinetic energy terms from the two-body, local, 
nonstiff potential terms has received considerable atten- 
tion [5lH53|, but we do not pursue this avenue here. 

In the context of nuclear motion only within the 
MCTDH method, it has been shown [27} that it is useful 
to decouple the orbital and A-vector equations for short 
times, which is enabled by the fact that the product of the 
inverse density matrix and the reduced operator p _1 W in 
Eq.® is in general slowly changing with time. Within 
the MCTDH implementation, p _1 W and the A-vector 
Hamiltonian H are taken as constant over a short time 
step, which approximation is denoted as "constant mean 
field" (CMF). The orbital and A-vector equations are 
integrated separately over the constant mean field time 
step, which is typically much bigger than, for instance, 
the time step used in the integration of the nonlinear 
equation for the orbitals. The error is determined by 
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backwards propagation and the CMF time step is ad- 
justed to keep it within a specified tolerance. 

We employ a similar method, but without intelligent 
error control and with constant stepsize, and with a pre- 
dictor/corrector scheme that appears numerically robust. 
The predictor step is identical to the CMF step in the 
MCTDH implementation, and the corrector step incor- 
porates a linear approximation to the matrices H(i) and 
/j -1 W(£). In the MCTDH terminology, this would be 
"linear mean held" (LMF). Thus, in the CMF predictor 
step starting at tg, with the matrices pq, W , and H for 
that time in hand, the wave function for the next time 
t\ = to + St is obtained as 



It = e- lHoSt A 



^c=(l-P(t)) 



h (D 



Pa 



'W 



t = tn...U 



(39) 



This first, predictor step in the propagation yields a first 
guess for A(ti) and c(t\), which yields first guesses for the 
matrices pi, Wi, and Hi at t\. These first guesses are 
used to propagate the wave function in a LMF corrector 
step, in which the first order Magnus approximation is 
used for the A- vector and a linear approximation for the 
product p _1 W is used. 

e -i(Ho+Hi)5i/2 j | Q 



Ax 
d _ 



+ 



(l-P(t)) 
t-t 



St 



c, 



' St 

t = to...t\ 



(40) 



The splitting of the orbital and A- vector equations over 
the mean field step is beneficial for, among other things, 
ensuring unitarity and parallelizing the algorithm. Al- 
though we have implemented versions of this method of 
higher order than linear (LMF), none exhibited its nearly 
unconditional stability. We use the expokit package [54J 
to calculate both the matrix exponential for the A- vector 
equation and the solution of the orbital equation. The 
exponential propagation of the orbital equation was the 
fastest method tried in this study, although we note that 
the explicit, basic Verlet method also gave good results. 

Our wave functions are Slater determinants and are not 
spin adapted; it is most efficient to calculate the high-spin 
case, so for a triplet we include projections of total spin 
M s = 1 , but therefore higher multiplcts are present in the 
configuration basis. However, we intermittently project 
the wave function on the proper spin subspacc to ensure 
that it is not contaminated by numerical error. A full 
description of the integration method will be presented 
in a forthcoming publication. 





Ro Nn £ elements 


Energy 


H 2 


1.4 9 14 3.0, 10.0, 10.0 
same with6> = 15° 1.1 xlO" 9 i 
same with6> = 30° 1.2xl0 -9 i 
HF limit 


-1.13362957146 
-1.133629573 
-1.133629572 
-1.1336295715 [55] 


Li2 


5.051 25 20 0.75, 3x 4.0 
eliptic basis HF 


-14.8715620178 
-14.8715619 [56J 


LiH 


3.015 21 19 1.0, 3x 5.0 
numerical HF 


-7.987352237 
-7.987352237 [57] 


CO 


2.132 21 19 1.5, 7.5, 7.5 
same withfl = 15° 1.1 xl0~ 8 i 
same with6> = 30° 8xl0" 8 i 
numerical HF 


-112.79090718 
-112.79090714 
-112.79090714 
-112.790907 [5JJ 


N 2 


2.068 21 19 1.5, 7.5, 7.5 
numerical HF 


-108.99382563 
-108.9938257 [57] 


N 2 


same basis, (14/10) CAS-SCF 
(14/10) Columbus cepvtz 
(14/10) Columbus ccpvqz 


-109.14184793(5) 

-109.132509251 

-109.140039408 



TABLE I: Converged Hartree-Fock energies from MCTDHF 
relaxation calculations and the FEM-DVR basis sets required 
to converge them, compared with literature values. For H 2 
and CO, the calculation is repeated with complex scaling with 
£o = 1, for two scaling angles. In these ECS results the last 
real digit and the imaginary components are not converged 
with respect to primitive basis. Also in the last entry, an 
MCTDHF relaxation calculation equivalent to a 10 orbital 
full CI MCSCF result for iV 2 , compared with results computed 
with cartesian Gaussian functions and a triple or quadruple 
zeta basis set. 



have been well established in the literature, others - in 
particular, the use of prolate spheroidal orbitals shared 
among all points in R — have not, and for this rea- 
son here we provide various ground, metastable, and 
excited vibrational state properties calculated with the 
present method. These are obtained by "improved re- 
laxation" (58l . l59j , in which the orbitals are propagated 
forward in imaginary time and the CI Hamiltonian (ex- 
cept when there is only one configuration) is diagonal- 
ized at every time step. We have also implemented a 
state-averaged version analagous to a state-averaged mul- 
ticonfiguration self-consistent field (MCSCF) calculation 
in which the orbitals are optimized to minimize the aver- 
age energy of the first TV cigenfunctions of the A-vector 
Hamiltonian. This procedure requires averaging the den- 
sity matrices and reduced operators for the first TV eigen- 
functions of the ^4-vector Hamiltonian and propagating 
their shared orbitals in imaginary time. 



VIII. GROUND ELECTRONIC STATES FROM 
IMAGINARY TIME PROPAGATION 



A. Fixed nuclei ground electronic states 



Of course, one requires initial state cigenfunctions to 
be used as a starting point for a time dependent cal- 
culation. While some aspects of the present method 



First, as a measure of the performance of the prim- 
itive basis in representing electronic wave functions, we 
list calculated fixed-nuclei Hartree-Fock energies for a va- 



Energy 



Natural occ. 



no. orbitals E (hartree) 


orb. 


B.O. 




N.B.O. 


1 


-0.5946128688 


1 


.99634 




0.99629 


2 


-0.5978526051 


2 


3.64 xlO" 


3 


3.68 xlO" 3 


3 


-0.5978974489 


3 


2.44 xlO" 


5 


2.48 xlO -5 


4 


-0.5978979622 


4 


2.36 xlO" 


7 


2.42 xlO -7 


5 


-0.5978979683 


5 


2.50 xlO" 


9 


2.6 xlO -9 


6 


-0.5978979683 


6 


2.94 xlO" 


11 


3.1 xlO -11 


exact 


-0.5978979686 


7 


4.4 xl0~ 13 


5 xlO -13 


Ref. [62] 


-0.5978979686 


8 


7 xl0~ 15 




1 xlO" 14 



TABLE II: Left: Energies of MCTDHF wave functions for 
ground state HD + as a function of orbitals, along with the 
exact result using our prolate spheroidal basis and the ex- 
act J = nonadiabatic result of Balint-Kurti et al [6^ . 
Right: natural prolate spheroidal occupation numbers for 
Born-Oppenheimer and non-Born-Oppenheimer HD + calcu- 
lations. 



riety of molecules in Table |TJ An MCTDHF calculation 
using full CI in the space of 10 orbitals space for N2 is also 
reported and compared with the corresponding calcula- 
tion using the Columbus quantum chemistry suite (60| 
and the correlation-consistent triple and quadruple-zeta 
bases of Dunning [6l[ . We also include results for a few 
molecules including complex scaling with £0 = 0, in which 



the coordinates of all electrons are continued into the 
complex plane. These latter complex scaling results were 
obtained by using the c-norm, not the hermitian norm, 
as explained in Section IVI1 and demonstrate that the 
electronic Hamiltonian has been accurately analytically 
continued to complex £. To achieve a given accuracy, 
these in general require slightly more DVR basis func- 
tions because of the oscillatory nature of the solutions 
under complex scaling. 



B. Nuclear motion: HD and natural orbitals for 
electrons and nuclei 

In our treatment the orbitals are used to span the 
entire range of internuclear distances R. Because the 
prolate spheroidal coordinates do not mimic the behav- 
ior of molecular orbitals - which asymptotically limit to 
atomic orbitals with constant size, whereas the prolate 
spheroidal coordinates continue to expand with increas- 
ing R - a greater number of orbitals is required to repre- 
sent a wave function with nuclear motion than one with- 
out. To precisely quantify this behavior, in Table [TTI we 
give ground state HD + energies calculated both with a 
numerically exact, converged calculation we performed 
using a large primitive FEM-DVR basis, and with the 
MCTDHF method for an increasing number of orbitals. 
One can see that sub-microhartree accuracy is achieved 
with three orbitals, and essentially the exact result is 
achieved with 5 orbitals. For HD + , our Hamiltonian is 
exact for J = and our exact result agrees with that of 





2 /C 






R 






W 


o)j 






FIG. 1: (Color online) Schmidt decompositon of the HD + 
ground state wave function: Natural prolate spheroidal (el- 
liptic) orbitals and conjugate natural orbitals in R for J=0 
ground state HD + , diagonalizing the reduced density ma- 
tricers in Eqs. (|4ip and (|42[) . corresponding to the occupa- 
tion numbers listed in Table [TTJ The origins of the prolate 
spheroidal coordinate system are denoted by black dots. 



Balint-Kurti [62| to all significant figures given. 

Also shown in Table [H] are the occupation numbers 
corresponding to the eigenfunctions of reduced density 
matrices for the ground J = state of HD + . These 
are the natural orbitals for electronic and nuclear motion 
in this coupled system. These natural orbitals and their 
eigenvalues provide a compact representation of the wave 
function known in the quantum information literature as 
a Schmidt decomposition (o8l - [70| . 



According to the theorems on which the Schmidt de- 
composition is based we may divide the HD + molecule 
with J = into two subsystems, namely that represented 
by (1) the coordinates of the electron and (2) by the nu- 
clear separation, R. Then the two reduced density ma- 
trices, that for electronic motion, 



Pel (e,v',<l>',Z,V,<l>)= / dRR 2 ^, V ^,R)r(e,v'^':R) 



(41) 
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V 


< E > 


(T) 


(V) 


in) 


<r?> 


(r2) 


irl) 


(R) 


(R 2 ) 


D 


5a Itt H 2 





-1.16008 


1.16006 


-2.32014 


1.5834 


3.1894 


1.5834 


3.1893 


1.4553 


2.1451 


0.0000 


8a Itt H 2 





-1.16088 


1.16085 


-2.32174 


1.5784 


3.1620 


1.5784 


3.1620 


1.4527 


2.1383 


0.0000 


Ref. [6J", [64J 6 


-1.16403"' 6 


1.16403 6 


2.32805 b 










1.4487" 


2.1270" 




5a Itt B.O. 









-2.32242 


1.5773 


3.1557 


1.5773 


3.1557 


1.4520 


2.1367 


0.0000 


5a Itt H 2 


1 


-1.14047 


1.14040 


-2.28087 


1.6270 


3.3691 


1.6270 


3.3691 


1.5482 


2.4787 


0.0000 


8<t Itt H 2 


1 


-1.14156 


1.14167 


-2.28324 


1.6295 


3.38286 


1.6295 


3.3826 


1.5482 


2.4817 


0.0000 


Ref. [63J 


1 


-1.14506 














1.5453 


2.4740 




5cr Itt HD 





-1.16164 


1.16158 


-2.32322 


1.57520 


3.14928 


1.57547 


3.15029 


1.44634 


2.11511 


-.0005391 


Ref. [65J 





-1.16547 






1.57119 


3.13009 


1.57148 


3.13120 


1.44223 


2.10432 




5o- lvr B.O. 









-2.32506 


1.5738 


3.1409 


1.5738 


3.1409 


1.4456 


2.1142 


0.0 


6-cr Itt LiH 





-8.03762 


8.03765 


-16.0753 


2.5808 


7.8354 


1.9864 


6.6936 


3.0834 


9.5398 


2.3458 


Ref. [66] 





-8.06644 






2.5651 


7.74517 


1.9719 


6.5857 


3.0610 


9.4197 




6-cr Itt BO 









-16.08629 


2.6219 


8.1238 


2.0066 


6.8593 


3.1285 


9.8404 


2.3306 



TABLE III: Properties of vibronic states. The H 2 calculation is from a state averaged calculation on the u — and v = 1 
states. Otherwise the energy of the ground vibrational state has been minimized. With six a and one ir orbital for LiH, fixed 
nuclei at 3.015, the dipole moment calculated was 2.2856 atomic units as may be compared with the prior result of 2.306 [67j . 



and for nuclear motion 

Pnuc(R, R) = 

/ d£dr)d(j) i/j(£,r),(j),R)ilj*(£,r),4),R') 

(42) 

have exactly the same same nonzero eigenvalues, pi. The 
complete wave function may be expressed in terms of the 
eigenfunctions, y>f(£, f], <f>) an d <Pi UC (R) of these matrices 
as 

1>(i, rj, <)>,R) = Y, Pl /2 vt& V, 4>) VT C (R) (43) 

i 

The pi are the natural occupations, and are a measure 
of the degree to which the parametric dependence of the 
prolate spheroidal coordinates upon the bond length fol- 
lows the change in the electronic wave function within the 
Franck-Condon region. In contrast, beyond pi the occu- 
pation numbers in cartesian coordinates x, y, z that do 
not follow the nuclei would be much higher. In Table |TT] 
we show two sets of occupation numbers, those for the 
Born Oppenheimer approximation to the ground vibra- 
tional state and those for the the numerically exact wave 
function whose energy agrees with Ref. (62[, and the oc- 
cupations are comparable. In Figure [T] we plot the pairs 
of corresponding natural orbitals in £, r\ (independent of 
cj> since m = 0) and in R, obtained from the exact wave 
function natural orbitals. In this example, only slight 
differences exist between these and those from the Born- 
Oppenheimer or improved adiabatic j43| wave functions 
in the same coordinate system. 

For a more complicated system, the concept of these 
coordinate-system-dependent natural occupations can be 
generalized. For a multiclectron wave function, the gen- 
eralization of the natural orbitals in R is straightforward. 
In this case, for the electronic degrees of freedom, we 



would have natural multielectron wave functions, not just 
orbitals, corresponding to the same set density matrix ei- 
genvalues. For a polyatomic system, we expect that the 
number of terms needed to converge the Schmidt decom- 
position analogous to Eq. (|43|) . as indicated by the the R- 
natural orbital or natural wave function occupation num- 
bers will be a measure of the suitability of a hypothetical 
geometry-dependent electronic coordinate system. One 
could compare two different choices of coordinate systems 
for the electronic degree of freedom (which like prolate 
spheroidal coordinates need not be orthogonal to R) by 
computing only the eigenvalues of the reduced density 
matrix, p nuc , for a suitable wave function. 



C. Vibrational states 

In Table IHlI we give properties calculated using MCT- 
DHF for the J=0 ground vibrational states of LiH, for 
H2, and HD, using a modest number of orbitals, with 
comparison to exact nonadiabatic results from the liter- 
ature. The reported values for H2 are obtained from a 
calculation in which the first two vibrational states are 
simultaneously optimized using the same set of orbitals, 
whereas the other results are from optimizing the ground 
v = state only. These calculations all use a single it or- 
bital and varying numbers of a orbitals. Differences in 
energies from the exact result on the order of several mil- 
lihartree for HD and H2 or tens of millihartree for LiH 
are apparent. The various expectation values differ by 
approximately one percent or less from their exact val- 
ues, even though in our multiconfiguration wave function 
a relatively small number of electronic orbitals have been 
used to span all gridpoints in R which number from 36 
for H2 to 48 for LiH. For the calculations Tables |TT] and 
IIIII and Fig.[2]we use nuclear masses run = 1836.152701, 
m D = 3670.483014, m Ll = 12789.395862. 
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FIG. 2: (Color online) Convergence of vibrational transition 
frequencies for Eb and LiH from state averaged MCTDHF cal- 
culations using one n orbital while optimizing 4 vibrational 
states for LiH and 10 vibrational states for H2. The differ- 
ence between the calculated transition frequency and the ex- 
perimental one is plotted with respect to the total number 
of orbitals, varying the number of a orbitals. For LiH, the 
errors using the Hartree-Fock Born-Oppenheimer curve are 
also plotted, with atomic masses as the arrows on the right 
side. For H2, the errors for the same transitions from a Born- 
Oppenheimer calculation with five <r, one 7r orbitals are also 
plotted as arrows to the right. 



In Tabic IIIII wc also report properties calculated for 
the Born-Oppenheimer wave function, i.e., the solu- 
tion xo{R)^Pi{r] R) obtained by diagonalizing the Born- 
Oppenheimer vibrational Hamiltonian using the ground 
Born-Oppenheimer electronic state ipi (r ; R) with orbitals 
optimized for each Ri separately, and using the atomic 
masses. For H2, the error in these results is comparable 
to that of our nonadiabatic MCTDHF wave functions. 
For LiH, we achieve better agreement with the previ- 
ously computed accurate values using the full nonadia- 
batic treatment than we do with the Born-Oppenheimer 
calculation using fixed masses. 

There is a striking observation to be made about these 
results concerning the convergence of our approach in 
which a single set of electronic orbitals is used for all R. 
By accounting for the nonadiabatic coupling terms in the 
Hamiltonian and using one small set of orbitals for all R, 
we achieve a better representation of the wave function 




0.0000836 




0.000182 




0.00446 





FIG. 3: (Color online) Electronic natural orbitals of ground 
state H2, with occupation numbers, from calculations with 
one 7T and six a orbitals. Left column: Born-Oppenheimer 
natural orbitals at R = 1.4 ao; right column: natural orbitals 
from nonadiabatic calculation of the v = state. 
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than wc do using the Born-Oppcnhcimcr wave function 
with orbitals optimized separately at each R. 

In Fig. [5] we show the performance of the method in 
representing the vibrational spectrum of H2 and LiH. For 
LiH, the first four vibrational states are simultaneously 
optimized using the same orbitals in these calculations, 
and for H2 the first 10 are optimized. The errors in the vi- 
brational transition frequencies are plotted with respect 
to the total number of orbitals. The corresponding er- 
rors in the transition frequencies for the vibrational states 
of the Born-Oppcnheimer curves (again, computed with 
atomic masses) are plotted as arrows on the right. These 
errors are comparable so the arrows overlap. 

Because in these calculations wc are simultaneously 
optimizing a set of vibrational states spanning a larger 
range of internuclear distances than v = and v = 1, 
the errors in Fig. [5] are greater than those in Table IIIII 
However, despite this fact the errors in the vibrational 
transitions may be made to be quite small even in a state 
averaged calculation. Wc obtain 4173.4cm -1 versus the 
correct value of 4161.1 for the v = to 1 transition of 
H a . 

In Fig. [3] we plot the natural orbitals from a Born- 
Oppenheimer calculation and the MCTDHF calculation 
with the same number of orbitals for the v = state, 
labeled by their occupations. Although some differences 
may be seen, particularly in the sixth a orbital, the over- 
all impression is that the two sets of natural orbitals are 
remarkably similar. That similarity suggests that the 
electron-nuclear correlation is substantially accounted for 
by the dependence of the orbitals on R via the prolate 
spheroidal coordinate system, since they have no other R 
dependence. A similar conclusion can be drawn from an 
examination of the natural orbitals for the state averaged 
calculation (not shown). 



IX. CALCULATION OF IONIZATION CROSS 
SECTIONS 



We calculation ionization probabilities and cross sec- 
tions using the flux formalism of Jackie and Meyer [7l[ ■ 
Three MCTHDF steps are involved: (1) relaxation to 
the ground initial state, (2) propagation from t = —T 
to t = in which a pulse of duration T is applied, and 
(3) propagation of \l/(0) forward in time until the ionized 
portion has been absorbed by complex the ECS grid in 
£. The wave function propagated during the third step, 
from t = onward, is saved and used in the following 
analysis. 

As per Ref. [7l|, the total ionized flux at energy E is 
defined as 



f(E) 



dt 



dt' (*(0) 



e i(H-E)t' P e -i{H-E)i 



*(0) 
(44) 



The flux operator F is defined as the flux through a hy- 
persurface, the region exterior to which corresponds to 
the breakup process of interest - in this case, ioniza- 
tion. Defining the heaviside operator Q(fi,f2, ■■■) to be 
the unit operator in this exterior region and zero within, 
the flux operator may be expressed as the commutator 
of the Hamiltonian with this heaviside function, 



F = i[H, 6] 



(45) 



Under appropriate assumptions and after some 
algebra [71|, one arrives at 



f(E) 



dt / dt' e 
Jo 



iE(t-t') 



(46) 

In this expression, the flux is obtained through matrix 
elements of the antihermitian part of the Hamiltonian 
between wave functions at different times; in deriving it, 
we have exploited the fact that the antihermitian part 
of the Hamiltonian lies only on the complex part of the 
ECS contour. 

In the present context, the antihermitian part of the 
Hamiltonian comes from the exterior complex scaling in 
the £ coordinates, and is nonzero only when at least one 
electron has reached the FEM-DVR element that termi- 
nates at £0, the start of the ECS tail. As long as £0 has 
been chosen sufficiently large, the antihermitian part is 
nonzero in the region corresponding to single or multiple 
ionization. In the MCTDH implementation for heavy 
particle motion, complex absorbing potentials (CAPs) 
are used instead of ECS in the exterior region to absorb 
the outgoing part of the wave function, and the validity 
of the above equations for the flux depend on the CAP 
being weak enough not to perturb the wave function in 
the inner region. In contrast, ECS is an analytic contin- 
uation of the Hamiltonian to complex coordinates, and 
does not perturb the solution in the inner region at all, 
unless a significant basis set error is present. 

In terms of the flux, the integral photoionization cross 
section (summed over final channels) is 



a(E) 



f(E) 



(47) 



where a is the fine structure constant, uj is the photon 
energy, u) = E — Eq where Eq is the ground state energy, 
and T(E) is the Fourier transform of the pulse from a 
length gauge calculation, for example, 



T(E) 



dt 8{t)e 



i(E-E )t 



(48) 



— T 



To evaluate Eq. (|4^|) it is necessary to evaluate the ma- 
trix element of H — W between wave functions at two 
different times, comprised of two different sets of orbitals, 
in an efficient manner. Although for the present appli- 
cations to H2 simpler methods would suffice, wc use an 
approach that will be applicable to larger systems. To 
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that end, we transform to a biorthogonal set of orbitals, 
in the spirit of the treatment of Malmqvist [72| . after 
which we may evaluate matrix elements of arbitrary op- 
erators, for instance the flux operator, via Slater's rules 
for zero, single and double excitations just as in the usual, 
orthogonal case. Thus, given orbitals and A-vectors at t 
and t' , we first transform the orbitals <p(t') into a new set 
ip(t') which obey a biorthonormality relationship to <f>(t): 
(<Pi(t')\</>j(t)) = 8ij. Whereas the MCTDHF orbitals 
are themselves orthonormal at all times, the ip functions 
alone obey no such relationship, 

j 

The full wave function at time t' has a new A-vector of 
configuration coefficients - which we denote as the B- 
vector, B - corresponding to its expansion in the new 
biorthogonal orbitals (p(f). In our notation the A-vector 
corresponds to the configuration basis \n(t')), and we de- 
note the configurations made from the <p(t') orbitals as 
|m), where (n(t)\m(t')) = 8^. We solve for B via 

#(*') =J2B. ffl (t')\n)(n\m) =^^(t')|n) 

mn n (50) 

A(t') = S(t')B(t') Snrn = (n(t')\m(t')) 

To solve these equations we must first construct Sft^, 
the overlap between configurations defined in terms of 
nonorthogonal sets of orbitals, a task which becomes in- 
creasingly more demanding as the number of electrons 
increases. We can take advantage of the remarkable fact 
that for full CI wave functions, although the matrix S 
is dense, its logarithm, In S, has sparse representations. 
In fact, the matrix InS is not unique, for the same rea- 
son that the multibranched complex function ln(z) is not 
unique, and it has both sparse and nonsparse represen- 
tations. 

The nonzero matrix elements of a sparse branch of In S 
may be calculated by applying Slater's rules for the ma- 
trix elements of a one-electron operator such as the ki- 
netic energy, which usually are applied to a determinant 
basis constructed of (bi) orthonormal orbitals. The ma- 
trix element In with respect to configurations \ n a (t)} 
and \rh(t')) is zero if these configurations differ by more 
than one index, and otherwise is given by Slater's other 
rules for a one electron operator as applied to any branch 
of the matrix logarithm Ins. 

For full CI wave functions, which are employed in the 
present work, the solution of the linear equation A = 
SB can be done in sparse arithmetic using the Krylov- 
space expokit [54| subroutine ZGEXPV, which performs a 
matrix exponential onto a vector. The solution is thereby 
expressed as 

B = exp(-lnS)A . (51) 

The transformation to the biorthogonal basis being 
done, we proceed to evaluate matrix elements of the anti- 
hermitian parts of the Hamiltonian operators appearing 
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FIG. 4: (Color online) Fixed-nuclei ionization cross section 
of H2, £ symmetry, calculated using the full wave function 
with one n and nine a orbitals, with a single pulse of inten- 
sity lxlCP 13 W cm -2 , frequency 1.1 hartree, and duration 
0.5 fs in the length gauge. Top: cross section from threshold 
(0.60449 hartree) to the 1 £„ threshold at 1.28 hartree. 
The thresholds are marked with arrows and the Fourier trans- 
form of the pulse is plotted with arbitrary units as dots. Also 
shown are the results of Sanchez and Martin [73|]. Bottom: 
magnification of resonance region, including result using only 
three sigma orbitals. Several other results are plotted that 
nearly coincide: lower intensity (1 xl0~ 15 W cm -2 ); velocity 
gauge; and more orbitals of sigma, pi, and delta symmetry. 

due to our use of exterior complex scaling, calculating 
orbital matrix elements and assembling them into the 
configuration matrix elements in the same manner as in 
constructing the A-vector Hamiltonian matrix H . 

The results for ionization of H2 in E symmetry (polar- 
ization parallel to the bond axis) are shown in Fig. @] We 
find that calculation is converged at a total of one tt and 
nine a orbitals. The other parameters of the calculation 
are given in the caption to Fig. 2J In obtaining Eq. 05] 
we assumed that (i/j(0)\i/j(t)} = (ip(0)\ijj(-t))* , and via 
backwards time propagation, we have verified that this 
identity is obeyed in general for these MCTDHF wave 
functions, at least to one part in 10 ~ 4 . To eliminate 
the oscillations of the Gibbs phenomenon in the Fourier 
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FIG. 5: Fixed-nuclei ionization cross section of Hb, £ sym- 
metry, using the treatment of Eq. ()52p to solve for only the 
perturbation to the wave function: results for one through 
four initial state a orbitals, with twice that number for the 
pulse and flux steps, using the same pulse parameters as in 
Fig- HI Circles: results of Sanchez and Martin [73| . 



transform over a finite interval we additionally multiply 
Eq.([46|) by a sinusoids in t and t', as cos |y , where T is 
the time for which we propagate the wave function af- 
ter the pulse. The result is converged to visual accuracy 
within approximately 1600 atomic time units. 

In Fig. 2] we also plot the results for one, two, and 
three a orbitals only. All reproduce the overall magni- 
tude and shape of the cross section, but are incorrect in 
the energy range where the autoionizing resonances ap- 
pear. The one-orbital treatment is featureless there, as 
the parent H^~ (£„) state is not represented in the basis, 
but otherwise correct. The two and three-orbital treat- 
ments reproduce the Fano lineshapes of the resonances, 
but place them at incorrect energies, and additionally 
their locations do not appear to converge until the nine 
a, one n result shown in black. 

However, the result converges to a cross section 
slightly different than the accurate results of Sanchez and 
Martin [73| . We may be able to understand this numer- 
ical behavior by realizing that at these intensities, only 
about l/1000th of the wave function has been ionized, 
and that these calculations have not converged that por- 
tion. This slow convergence behavior would seem to be a 
problem for the utility of the MCTDHF method for de- 
scribing photoionization or other perturbativc processes 
in general. One would like to treat perturbative problems 
just as well as nonperturbative ones, but the variational 
ansatz of the MCTDHF wave function will use the varia- 
tional flexibility in the calculation to optimize the larger, 
unperturbed (initial state) portion of the wave function 
at the expense of the smaller components in which we 
are more interested. 

In the limit of a large number of orbitals, the MCT- 
DHF equation should converge to the exact result. It 
is likely that this number is much larger than we have 



used in the calculations shown in Fig. [4] as additional 
orbitals are likely to mostly further optimize the corre- 
lation within the ground state until enough have been 
added so that the occupation numbers of the natural or- 
bitals describing ground state correlation have fallen at 
least two orders of magnitude. There is, however, an 
alternative approach. 

X. MORE EFFICIENT CALCULATION OF 
IONIZATION 

The problem that additional orbitals mostly serve to 
improve the description of the initial state is particular to 
the present application to calculate a perturbative result, 
and more intense-field applications would not suffer from 
it. This state of affairs in unsatisfactory, but fortunately 
there is a straightforward solution to this problem. We 
can calculate not ^(t), the full wave function, but the 
quantity we label ^'(t), the change in the wave function 
due to the pulse: 

#(£) = e- ?;Bot ^(0) + $'(t) 

d (52) 
i— = H{t)$>'(t) + V(t)e- tEot V{0) 

where Eq is the initial state eigenvalue and V(t) is the 
perturbation. This formulation modifies the MCTDHF 
working equations, Eqs.® and (|13p . introducing driving 
terms to each, but we were not immediately able to im- 
plement the orbital driving term in a numerically stable 
way. We thus double the orbital dimension at the start 
of the pulse, generating additional orbitals by operating 
with the dipole operator upon the occupied orbitals of 
the initial state, and causing the orbital driving term to 
equal zero at the start of the pulse. We then calculate 
ty°(t) = e iEat ^'{t) by modifying Eq.([T3]) as 

i^-A°it) = (H(t) - E )A°(t) + V(t)A(0) , (53) 
ot 

where the matrix V is the matrix of the dipole opera- 
tor in the nonorthogonal basis of orbitals at times t and 
0, V )1r -, = (il(t)\fi\fi'(0)). This equation is solved with 
routine ZGPHIV in the expokit package [54) . 

The results are significantly better than in our first 
treatment, calculating the entire wave function, although 
we find that we cannot accurately calculate the entire 
range between the H2 ionization thresholds with the sin- 
gle pulse we used for the calculation of the full wave func- 
tion Vf'(t). We focus on the resonance region, where our 
hui = 1.1 hartree, 0.5 fs pulse is centered. The results 
are again insensitive to intensity across the whole energy 
range. 

In Fig.[5]we show that the ionization cross section con- 
verges to essentially its correct value using a orbitals only. 
In our treatment we double the number of orbitals going 
from the ground initial state to the propagation steps 
of the calculation; thus the minimum is two. We can 
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see that the minimum one orbital (Hartrce-Fock) ground 
state, two propagation orbital treatment yields an ioniza- 
tion cross section without the correct resonance features; 
the two orbitals of the propagation of ^'(t) correspond 
to the ground H^~ la g cation state and the wavepacket 
ionized in its field. The resonances, which are based on 
the la u cation state, are thus not represented, initial, 
six propagation orbitals, the resonances appear in essen- 
tially their correct locations. Two additional propagation 
(one additional initial) orbitals are all that is needed to 
give good agreement with the calculations of Sanchez and 
Martin [z3|. 

XI. CONCLUSION AND OUTLOOK 



FEM-DVR, see Ref. jUj]. We evaluate all the integrals 
by quadrature; the only non-polynomial terms approx- 
imated by quadrature are the repulsive inverse integer 
powers that appear, for example, in the centrifugal po- 
tentials. 

We refer to the matrix elements using the notation of 
Eq. ([2"5)) representing the Hamiltonian we use here, which 
is exact for J = except for the omission of the two- 
electron terms in P. The Hamiltonian is otherwise ac- 
curate except for Coriolis coupling for J / 0. In the 
equations below, f££ refers to an electronic FEM-DVR 
product basis function of Eq.(|19p. \ to a primitive DVR 
function of Eq. ([T5)) . and x' to its first derivative. 

The matrix elements in R are straightforward, 



We have explored the formulation of the MCTDHF 
approach both for fixed nuclei and including nuclear mo- 
tion for application to any diatomic molecule within the 
nonrclativistic approximation. Furthermore methods for 
overcoming several important technical barriers to such 
calculations have been demonstrated. The use of pro- 
late spheroidal coordinates, and an expansion of the wave 
function including nuclear motion in terms of configura- 
tions of orbitals which depend on R only through the 
dependence of the underlying coordinates on internu- 
clear distance are crucial parts of the strategy described 
here. We have demonstrated that the use of such or- 
bitals gives a rapidly convergent representation of the 
vibrational states of diatomics in calculations that avoid 
the Born-Oppenheimer approximation. Furthermore, we 
have shown that photoionization cross sections can be 
extracted from the MCTDHF wave functions in full di- 
mensionality, and demonstrated that accurate photoion- 
ization cross sections can be calculated using a method 
for solving for only the perturbation caused by a time- 
dependent potential. The methods we describe here 
should be immediately applicable to calculation of the 
single (or total) ionization of larger molecules, as well as 
studies of Auger spectra with or without nuclear motion. 
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Appendix A: One electron matrix elements 

We follow Ref. [33| in the derivation of the matrix 
elements of the electronic Hamiltonian in our prolate 
spheroidal DVR basis, and provide additional formulas 
for the nonadiabatic terms. For details on the use of 



where we use the shorthand p = \J (£ 2 — 1)(1 — rf-\ may 
be expressed using the identities 




16 



-iei 



£ - 1 + (Cc + q^) 2 (C 2 - 1) 2 + p 2 (v + O 2 



2/i e 2fi R 



i-4 + (vk + < a f (i - vl) 2 + p 2 (£ + an) 2 



action of the Jacobian of Eq. (J8j) onto the orbitals within 
the mean field step. 

We begin by defining regular and irregular Legendrc 
functions [73 with an additional normalization factor, 



2/i e 2fj, R 
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so that the Neumann expansion of I/V12 may be written 
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Thus to compute the two-electron integrals we require 
the matrix elements of Pi m (£,<)Qim(£,>) in the DVR ba- 
sis in £. This function is the Green's function for the 
following equation: 



c 2 



In the present applications the pulse was parallel to the 
bond axis, and the dipole operator is given in length and 
velocity gauge as 
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respectively. 



Appendix B: Two-electron matrix elements 

We can follow the method of Refs. [HI, HH to con- 
struct the matrix elements of l/r 12 in the DVR basis. 
We evaluate the matrix elements of a multipole expan- 
sion of this operator by solving Poisson's equation in the 
radial (£) coordinate while evaluating the rj integrals by 
quadrature. The resulting matrix elements retain the 
sparsity of the DVR representation, being diagonal in £ 
and r\ and off-diagonal only in the M. quantum numbers 
of the electrons. This point is crucial for the present 
time-dependent application, as it means that the two- 
electron transformations do not dominate the computa- 
tional time, which is instead primarily determined by the 



where W is the Wronskian of the two Legendrc functions, 
which has the value 

(-l) m 2 2m r ( l + m + 2 ) p ( i+m+i \ 

Expressing the operator in square brackets in Eq. (|B3p in 
the DVR basis and approximating the matrix elements of 
both sides of that equation using the DVR quadrature, 
we arrive at an expression for the matrix elements of 

fim(£<)Qjm(£>)> 
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(B5) 

In Ea. (|B5|) . N is the last Gauss- Radau gridpoint in £, 
corresponding to a DVR function discarded to enforce 
the correct boundary condition at the end of the £ grid 
on the solution of Eq. (|B3|) . and T; m is the the ma- 
trix of the operator in square brackets in that equation 
in terms of the DVR functions Xi(£) (f° r even m) or 
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\/(£ 2 1)/(C 2 ~~ l)Xi(0 (f° r °dd m), which are normal- 
ized to one with respect to integration over £. Those 
matrix elements are 



(Tlm)ij — — Si 



£-1 



+ /(/ + !) 



(B6) 



Using the expression for R 1 ^ in Eq. (|B5I) we obtain 
the final result for the two electron matrix elements in 
our DVR basis, 



which has exactly the form in Eq. ([20|) . This expression 
depends also on having used a fixed DVR quadrature 
to approximate the r\\ and 772 integrations, and a given 
quadrature order cannot be used for arbitrarily large I 
values in the Neumann expansion in Eq. ([B2[) . In our 
numerical calculations we use l max = N v where N v is 
the number of DVR functions in rj. We have found 
this choice to be optimal when using this algorithm 
implemented for spherical polar coordinates, and at least 
near optimal for the present case of prolate spheroidal 
coordinates. 
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